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Array Calibration in the Presence of Multipath 

Amir Leshem, Member. IEEE, and MaU Wax, Fellow. IEEE 



Abstract — We present an algorithm for the calibration of sensor 
arrays in the presence of multipath. The algorithm is based on 
two sets of calibration data obtained from two angularly separated 
transmitting points. We show the similarity between the caUbra- 
tion problem and blind identification of SIMO systems and analyze 
the identifiability of the problem. Simulation results demonstrating 
the performance of the algorithm are included. 

index Terms—Arrtty calibration, blind channel identification, 
DOA estimation, multipath. 

I. INTRODUCTION 

MODERN super-resolution direction finding techniques 
such as minimum variance [I], MUSIC [4], subspace 
fitting methods [8], and maximum likelihood [12] presume the 
knowledge of the array response. 

As the analysis of these techniques show [6], [7], [ 1 1 ]. any in- 
accuracy in the presumed array response results in severe degra- 
dation of performance. The measurement of the array response, 
which is referred to as array calibration, is therefore a crucial 
step in the implementation of these techniques. 

The existing calibration techniques [5], [9] are based on mod- 
eling the array response by a free-space model perturbed by 
an unknown coupling matrix and sensor location uncertainty. 
These unknown pajameiers are estimated together with the un- 
known signal parameters, assuming known or unknown source 
location. Yet, for general arrays with arbitrary sensor responses, 
these methods are no longer adequate since these modeling as- 
sumptions are no longer valid. 

In this paper, we address the problem of measuring the array 
response of arrays with arbitrary sensor response in the pres- 
ence of multipath. This problem is important since multipath is 
essentially unavoidable in practice, and it sets the limit on the 
achievable calibration accuracy. 

The organization of the paper is as follows. In Section II, we 
formulate the problem. In Section III, we present the proposed 
solution. In Section VI, we present simulation results demon- 
strating the performance of the algorithm. Finally, in Section 
VII. we present some concluding remarks. In Section V, we con- 
sider the similarity and the differences between the calibration 
problem as presented here as well as the problem of blind iden- 
tification of multiple FIR channels. 
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Fig. I . Calibration setup (with one reflectaon)< 

II. Problem Formulation 

Let a{$) denote the p x 1 vector of the array response to 
a source impinging from direction 0. The array calibration 
problem amounts to measuring a{0) for 0 € (0, 27r). It is 
usually performed by transmitting a signal from some location, 
rotating the array, and measuring the array response at each 
angle. Unfortunately, in many cases, the measured response 
is composed not only of the direct path from the transmitting 
point but also of multiple reflections from the surroundings; 
see Fig. L In the case of arbitrary array response, we can no 
longer resolve the multipath from a single set of measurements 
since the measured data can be considered to be the "true" 
array manifold. This situation is similar to the problem of blind 
identification of SIMO systems, wherein without any a priori 
knowledge of the signal, a single channel is not identifiable, 
and two channels are identifiable, even using second-order 
statistics only. 

To cope with the multipath problem, we propose to carry 
out the calibration twice, i.e., rotate the array and measure 
the received array vector as a function of 0^ yet each time 
use a diflferent transmitting point. Let y,(^) denote the p x 1 
vector received at the angle 0 from the lih transmitting point 
(/ = 1; 2). Assuming that the reflections are considered as 
point sources and all multipath effects are completely coherent 
with calibrating signals, i.e.. each path diflers by a complex 
reflection coefHcient from the direct path, we get 

i=I,2 (1) 

where 

I direction of the tth reflection in the lih set; 
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Pi^ I complex coelTicient rq)resenting the phase shift and 
the amplitude of the ith reflection in the /th set; 

ri number of reflections in the /th set; 

ni{0) noise vector for the angle 0 in the /th set. 
Since the array manifold is measured relative to some arbitrary 
point and the relative angle between the measurement points is 
known, we can assume without loss of generality that p\^\ = 1 
and ^1, 1 = ^1, 2 = 0**. In addition, since the reflecting objects 
remain flxed while the transmitting point change, the relative 
directions of the reflections are different, i.e., i # 2 (i / 
1). 

Assuming that the calibration process consists of N measure- 
ments taken uniformly on ^ € [0, 2ir), it follows from (1) that 
the measured data is given by 

..(t)=|:«."(^*«''0*"'« 

/ = l,2;A: = 0,....yV-l (2) 

where we use n/(fc) to emphasize that the noise is not angle 
dependent. Note thai we have included the direct path with the 
multipath. 

In our solution, we make the following assumption: 
Al) All reflections Oij are a multiple of the basic rotation 

27^/^^ 

Assumption A I ) serves as a very good approximation when the 
grid is fine. 

The array calibration problem can now be formulated as fol- 
lows. Given the two measured data sets 



estimate the array manifold 



III. THE Maximum Likelihood Estimator 
The proposed solution is based on two steps: 

i) estimating the reflections* parameters {pj, Bi)\ 

ii) estimating the array manifold using the estimated reflec- 
tions; 

where pi = [/>!.«. ■ » />n,/] is the vector of the re- 
flection cocfTicicnts at the ^th set of measurements, and 
Oi = (^11, . ■ , By^j] is the vector of the reflections* DOA's 
at the /th set of measurements. 

For the first step, we have two approaches. The first approach 
uses the LS estimator, which is identical to the MLE under the 
assumption of while Gaussian noise. This estimator is derived in 
this section. The second approach uses a simplified LS, which 
we derive in the next section. 

The second step is derived by a least squares solution, which 
is the MLE under the assumption of white Gaussian noise. This 
step is performed identically in the two approaches using the 
results of the first step. 
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To cany out the derivation of the MLE, let wi denote the TNA x 1 
vector whose A:th element is p*, i if 27f k/N = i and zero 
otherwise. Mathematically, this is expressed as 



(3) 



where 6{$) is the deha function. Let Om be the A' x 1 array 
manifold of the rath sensor 

»„<0).*.(|),-.^(?^^)]'<4, 



Om = 



(5) 
(6) 



and VfftjW and nju,i{k) are the mth element of yi(27r/:/iV) 
and Hi (A:), respectively. 
With this notation, we can rewrite (2) as 

n 



53«;,(A:i.,)n... 

>«0 



On* +n/ 



(7) 



where kij = {Bijl{2r /N)), and Pk is a permutation matrix 
that rotates the zeroth element of a into the /cth position defined 
by 



- / ^» ifm + A: = norm+ - TV = n 
(/'fc^m.n - otherwise. 



(8) 



The last equality in (7) is due to the fact the w{k) = 0 if /3i 
such that 2iTk/N = ^<,i. Denoting 



(9) 



it follows that Ai is mN xN circulant matrix generated by wi 



A, = 



••• w,{N-2) w,{N-l) 
wi(N-l) i«,(0) • • v„{N-2) 



(10) 



••• w,iN-l) «;,(0) 

Thus, we can rewrite (2) as 

V„ , = AjOm+nm.i l<rn<p. (II) 

Since At is a circulant matrix, it is diagonalizcd by the DFT 
matrix of order N, and its eigenvalues are given by the DFT of 
the generating vector ti/j [2]. Therefore 



F"A,F = diag{Fii»i> = diagfwi} 



(12) 
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where F is the normalized DFT matrix of order A' {FF^ = /), 
and \oi = Fwi is the DFT of wi given by 

""(^) = ^ E PM'^-^*^'''''/''^"'''- (13) 



which can also be rewritten as 



I) 



Hence 



>l, = Fdiag{w,}F". (14) 



With this representation of the matrices Ai^ we can derive a 
somewhat simplified expression of the MLE. 
Ut 



I 

ym,2 



and A 



Assuming that the noise is white and Gaussian from ( 1 1), the 
MLE is given by 

|d|, d,„ ^i, 6*2) Pn Pj] 

p 

= ajg inin V liy"*-^"!!'^- 



Minimizing first with respect to a„», we obtain 

ii,„ = {A»A)-'A'*v„. 
Now, from the defmition of A and (14), we obtain 



(16) 



Notice that this estimator involves all the reflections parameters, 
i.e., the DOA's and the reflection coefficients, in a highly non- 
linear fashion and, hence, is computationally unattzactive. 

rv. Simple LS Estimator 

In this section, we derive a simplified LS estimator for the re- 
flections DOA*s and reflection coefficients. This estimator, to- 
gether with the estimator for Om given in (18), consists of the 
simplified LS estimator for the array manifold. 

Substituting (14) into ( 1 1 ), we obtain 

F^on. = diag(ti;,)-Uyn».i " f'hn.t) (22) 

where hmj = F^^n^^i. Since this holds for both sets of mea- 
surements, we obtain 

W^rnVm.l -«m.l) = '(ym.2 -^.2) (23) 

which can be rewritten as 

°ym,2 " ^2 Oym.l = ^1 ^^rn,2 - ^2 ^ »Wn, i (24) 

where o denotes eiementwise multiplication. 

Since the right-hand side of (24) is "noise," a possible LS 
estimator for the reflections' parameters is given by 



(17) 



Substituting (13) into (25) yields 
[^1, ^2, Pi, P2] 



(25) 



where Wi = diag(Ffi;i). Substituting (14) and (17) into (16) 
yields 



where 



Jb-0 



(19) 



(26) 



and y„, / = F"y„^ t. Finally, substituting (14) and (18) into 
(15), we obtain 



[^1, ^2, Pi'. P2] 
= arg ^ inin L,„,, - FW^DWl^F^y^X 



r 2 



(20) 



where/:*, I = {ei^i/{2vk/N)). Denoting 

Ti(Ar,6^)r=y^.i(^)c--'** (27) 

T2(fc, 9) =i/n..2(A:)c-'*' (28) 

and (29), shown at the bottom of the page, we can rewrite (26) 
as a linear problem in p (recall that we have assumed pi, 1 =1 
and^i.i = 0**) 

f^i, p] = inin ||B„.(ai, ^2)P -H y.„.2ir'- (30) 



T2(0,<92.i) • • Y2(0.^n.l) -Ti(0,i9i,2) "Ti(0,tf^,.2) \ 

,T2(/yr 1,^2.1) T2(yV - !,<?,...,) -T,(.V~1,<9,,2) -TiiN - i, 9,^,2) J 

P = (P2,l. ■•1 Pn, ii Pl.2i •••» Pra.'i)^ 



(29) 
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This estimator is based on the data of the mth sensor only. 
Clearly, we can improve the performance by combining the in- 
formation from all sensors. This yields 

p 

f^i.^2,pl =arg imn V ||fin.(ei. »2)p + ym.2ll^- 

(31) 

To evaluate this estimator, we first rewrite it in matrix form as 

(32) 

(33) 



\0i, flj, p] = arg ^iiun \\B{Bu *2)p + y2ll^ 



where 



and 



BiOu $2) = [BliBu «2), ■ • • . BliOu 62)]' 



y2 = [vn^ " ' ^ yp2] ' 



(34) 



Minimizing first with respect to p, with $1^ B2 being fixed, we 
obtain the well-known least squares solution 

p = - (B(eu e2)"B{Bu B2)Y' B{Ou 02fif2' (35) 
Substituting (35) back into (32), the resulting estimator of the 
direclions-of-arrivai of the reflections is 

[^i. ^2] = arg imn \\PD(9uB.)iym,2)f (36) 

where is the projection on the orthogonal comple-- 

menl of the subspace spanned by the columns of B{9i , ^2) 

pX 

= / - B{$u Bi) {B^'iOu e2)B{eu tf2))"' b"{9u e^y 

(37) 

The structure of this estimator is similar to that of the deter- 
ministic signal maximum likelihood DOA estimator. Hence, the 
optimization methods developed for this problem, including the 
alternating projections [12] and the clustering methods [3] can 
be used. 

With the estimated parameters at hand, we can use (16) to 

estimate the array manifold. First, we obtain an estimate w of tii 
by substituting the /j's and the ^'s into ( i 3). We then get 

By substituting W - diag(wi) into (14), we obtain 



which when substituted into (16) yields 



(39) 



where 



(40) 



(41) 
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V. The Relation to the SIMO Blind Equalization and 

iDENTIFIABnJTY RESULTS 

In this section, we cast the calibration problem as the identi- 
fication of a single input multiple output (SIMO) system. This 
will enable us to derive identifiability conditions, as well as 
present an alternative derivation of the LS estimator for reflec- 
tions parameters. 

First, note that we can rewrite ( 1 1 ) as 

V^.i^^i^'O^'^'^.i ' = 1.2 (42) 

where 



lo, 



if^t.< = 



27r(A^ - fc) 



N 



(43) 



otherwise. 



That is, the measurements are just a spatially filtered version of 
the "signal" by FIR fikers with coefficients pi, x at Bi^i and 
zeros otherwise. Our problem can now be stated as follows. 

Given the output of two linear systems driven by the same 
signal reconstruct the input signal. 

The problem is in the form of blind identification. However, 
several differences between our problem and the conventional 
blind identification problem exist. 

1 ) The signal is periodic with known period. 

2) The measurements are taken along a single period. 

3) We have several pairs of output signals: one for each ele- 
ment of the array. 

4) The filters are sparse, i.e., most of the coefficients are 
zero. 

5) The length of the filters may be the same as the number 
of samples. 

We next develop the LS estimator as a natural variation on the 
method of [10] in the frequency domain. 

Using the convolution theorem (remembering that our signal 
is periodic), we obtain 

where o denotes element-wise multiplication. Hence 

ym.i =ym,2«>^i- (^5) 

After some^lgebraic manipulations, using the relation between 
ti; and h, we obtain the noiseless version of (24). 

The fact that our LS estimator can be derived using the ap- 
proach of [ 10] enables us to give a sufficient condition for iden- 
tifiability. This condition is obtained by translating the sufficient 
condition for identifiability of [ 1 0]. However, since our channel 
is sparse, we will be able to obtain stronger identifiability con- 
ditions. 

To that end, note that the problem is identifiable for channels 
with signature (ri, T2) (i.e., has a unique solution with at most 
Vi reflections in the first set of measurements, and, at most, r2 
reflections in the second set in the noiseless case) if (32) has a 
unique solution with the true (^1 . ^2) while having no solution 
with any other substitution of tf'i , ^. For this condition to hold, 
it is suflicient and necessary that for any pair {0\, $2) where 
^1 € C^" and ^2 € C^''^ the matrix , ^2) has full column 
rank. We shall elaborate on this to obtain some further condi- 
tions, which will be easier to verify. To simplify notation, we 
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will work with a single matrix Bk instead of with the full ma- 
trix B, The generalization is straightforward though notation- 
ally complex. 
Ut 



ym.l(0) 



(46) 



and (47). shown at the bottom of the next page. Note that 

Bmi0i,02) = [Ai\7]C{eu02y (48) 

The second matrix is always full column rank due to the Van- 
dermonde structure of each block. Thus, the identifiability con- 
dition boils down to having the fu^t matrix preserve the column 
rank. Similar to the condition in [10], we can now split this 
condition into two conditions. The first demanding informative 
array manifold, and the second is a condition on identifiable 
channels. Factoring A/ similarly to (7), we obtain 



A, = 



(49) 



Thus, the following immediately follows. 

Theorem 5.1: Let L = 2{ri -f r2 - 1). Assume that for 
every 0 < fc < I. - 1. ttni(^) 0 and iiJi(^) and W2{k) are 
not simultaneously zero; then, there is a unique solution to the 

problem (30), 

Note that our conditions depend on the number of reflections 
rather than the channel length. We can, of course, weaken the 
condition above. However, the above condition for the array 
manifold typically holds, leaving us with conditions on the 
channels that hold, e.g., if the channel polynomials do not have 
conrunon zeros on the unit circle. 

VI, Simulation Results 

In this section, we present the results of several simulated 
experiments that demonstrate the performance of the algorithm. 

In all experiments, the array consisted of two sensors that were 
2.5 A apart . and the number of reflections was 2, i.e., = = 
2. 

In the first experiment, the directions of the signals 
were 9i,i = 61^2 - O"", and those of the multipath were 
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Fig. 2. Anay manifold errors versus SNR. Multipath conditions: = 1 , 
pi a = 0.03 + 0.05i,/>3., = l,/»a.a = 0.13 + O.I97. , = = U^ 
013 = IS**. B22 = 40*. S/Mi = 25 dB, S/M^ = 13 dB. Solid Unc: error 
after application of (he algorithm. Dashed line: error due to multipath (first set 
of measurements). 



^1,2 = 15**, 02,2 = 40** (Note that this does not limit the 
generality of the simulations since, as explained earlier, we 
can align the direct paths of the two measurements and only 
estimate the angles of the multipaths relative to the direct 
path. Typically, after alignment, the multipath will arrive 
with different AOA*s, due to the fixed geometry of the re- 
flectors.) The reflection coefficients were p, i = 2 = 1, 
pi^2 — 0.03 + 0.06i, P2.2 = 0.13 -h 0.19j, which corresponds 
to signal to multipath ratios of 25 and 13 dB, respectively. At 
each SNR. we have performed 25 Monte Carlo trials. Fig. 2 
shows the array manifold error averaged over all DOA's as a 
function of the SNR. While the solid line represents the error 
after application of the algorithm, the dashed line presents the 
array manifold error in the first set of measurements. The array 
manifold estimation MSE is computed by 



MSE = 



1 

NM 



where M is the number of trials, and is the grid size. In all 
experiments, Af = 25. AT = 360, and a(27rA;/A^, t/i) is the ntth 
estimate of a(27r/:/yV). To gain further insight into the perfor- 
mance of the algoritlun, we present in Fig. 3 the results of a 





0 










0 





(47) 
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Fig, 3. Airay manifold errors versus DOA, Multipath conditions; (h,t = 1 , 
Pi a = 0.03 + n.05j.pa.i = l.^a.a = 0J3 + H.lOj. /?u = l^ai = 0**. 
^>t, = ^jj = 40 ». S/Afi = 24 dB, S/Afa = 13 dB. Bottom line: error 
after application of the algorithm. Upper lines: Error due to multipath (in two 
sets of measurements). 
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Fig. 4. Array manifold errors versus SNR. Multipath conditions: pi, i = I, 
Pi a = 0.21 + 0.05>.pa.i = I.Pa.a = 0.13 + 0.19>.dn = B^i = O*'. 
Bi-j = 15'^. $23 = 2.j'. S/Mi - 13 dB. S/Afa = 12 dB. SoUd line: error 
after application of the algorithm. Dashed line: error due to multipath (first set 
of measurements). 



single experiment performed at SNR of 50 dB. We can see that 
the error aAer (he application of the algorithm is much smaller 
than the error at the raw set of measurements. Moreover we see 
that the error is about the same for all DOA's. 

In the second experiment, the directions of the signals 
were = 0i^2 - O"*, and those of the multipath were 
^^.2 " ^2.2 - 25**. The reflection coefficients were 
p^^ = ^1.2 = h Pi. 2 = 0.21 +0.05j\ P2,2 = 0.13 + 0.19; 
corresponding to signal to multipath ratios of 13 and 12 dB. 
At each SNR, we have performed 25 trials. Fig. 4 shows the 
array manifold error averaged over all DOA*s as a function 
of the SNR. Whereas the solid line represents the error after 




Fig. S. Array manifold enors versus DOA. Multipath conditions: pi . i = 1 , 
p,.a = 0.21 f 0.05j\p,,i = l.pa,2 = 0.13 + ClOj. tf,, = fiji = 

= 15% Bm = 28«. 5/M, = 13 dB. S/M^ = 12 dB. Bottom line: error 
after application of the algorilhro. Upper lines: eirar due to multipath (in two 
seta of measurements). 




Fig. 6. Error as a ftmction of ^aa - A^ia. Multipath conditions: p, . i = I . 
pi a = O.U+0.2>,pa.i = 0.7, pa. a = 0.055 +0.2>. SNR = 30 dB. Dashed 
line: the eiror due to the multipath. Solid line: the error after the application of 
the algorithm. 



application of the algorithm, the dashed line presents the array 
manifold error in the first set of measurements. Fig. 5 shows 
the results of a single experiment performed at SNR of 50 dB. 
We can see that the error after the application of the algorithm 
is much smaller than the error at the raw set of measurements. 
Moreover, we see that the error is about the same for all DOA's. 

In these two experiments, we clearly see that the improve- 
ment is not only obvious, but the error reduces to the level of 
the measurement noise. This demonstrates that the multipath is 
completely removed. 

In the third experiment, the relative angular separation 
between the reflections in the first set of measurements was 
held fixed at i 5**, whereas the relative angular separation in the 
second set of measurements varied from 19-51'' in steps of 4". 
The reflection coefficients were pi, i = I, ^i, 2 = 0. 1 1 + 0, 1>, 
P2,i = 0.7. and p2,2 = 0.055 + 0.2; , and the SNR was 
30 dB. The results are presented in Fig. 6. Notice that the 
performance of the algorithm is essentially independent of the 
angular separation. 

VII. Concluding Remarks 

We have presented a novel method for the calibration of 
sensor arrays in the presence of multipath. The method is based 
on measuring the array manifold from two angularly separated 
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locations and involves a solution, of a multidimensional opti- 
mization. The method does not depend on the relative angular 
locations of the reflections. 
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